//public head 
#ifndef XWB_MATH_FIT  
#define XWB_MATH_FIT  

#include <opencv2/core/core.hpp>
#include <opencv2/imgproc/imgproc.hpp>
#include <opencv2/highgui/highgui.hpp>

#include <vector>  
#include <cmath>
#include <iostream>

#include <time.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>

typedef struct mt
{
	double *data;
	int m;
	int n;
}matrix;
#define STH 1// 0 meanpk, 1 bgk, 2 ost, 3 threshold value
const static float noisecounter = 0.4f; //connected region threshold value
const static float snrth = 1.0f; // snr threshold value

const static float gsegma = 0.9f;// gaussian kernel


#ifndef PI
#define PI 3.1415926535898
#endif

#define FFT3 //use fftw3

#define JGGL 0.2 // dye threshold

#define RPIXEL 1.05 // register precision

#define TPS 2.05 // template size

#define TPLC 4 // template cycle

//#define CHT // move template 

#define FBI 9

#define OPMP // is open open_mp

//#define SUBFFT // fft registration accuracy

#define PFT 0.4 // good peak threshold

#define SMV 10.0 // score map parameters

#define STV 4.0 //good score threshold

namespace xwb{
	////////////////////////////////////////////////////////////////
#define EXIT_SUCCESS 0
#define EXIT_FAIL 1 

#define EL(A,i,j) A.data[i*A.n+j] // a way to access element i,j of a matrix. In the above example, EL(A,1,0) will access A[4]
	// A macro to allocate memory for a matrix and checking if the allocation was succesful
#define NEW_MAT(A) A.data = (double*)calloc(A.m*A.n,sizeof(double)); if (A.data==0) {/*printf("MEMORY ERROR");*/ return EXIT_FAIL;}
	// A macro to see if a function return a failure. Handy to implement a cascading failure
#define EVAL(F) if (F==EXIT_FAIL) {return EXIT_FAIL;}


	static void printMat(matrix A)
	{
		int i, j;
		double toPrint;

		for (i = 0; i<A.m; i++)
		{
			for (j = 0; j<A.n; j++)
			{
				toPrint = EL(A, i, j);

#ifdef OFFLINE
				printf("%.4le ", toPrint);
#endif

			}

#ifdef OFFLINE
			printf("\n");
#endif

		}

	}

	static void fprintMat(FILE *fh, matrix A)
	{
		int i, j;
		double toPrint;

		for (i = 0; i<A.m; i++)
		{
			for (j = 0; j<A.n; j++)
			{
				toPrint = EL(A, i, j);

#ifdef OFFLINE
				fprintf(fh, "%.4le ", toPrint);
#endif

			}

#ifdef OFFLINE
			fprintf(fh, "\n");
#endif

		}

	}

	static double vectorMagnitude(matrix A)
	{
		double sqrSum = 0, ai;
		int i;

		for (i = 0; i<A.m; i++)
		{
			ai = EL(A, i, 0);
			sqrSum += pow(ai, 2);
		}

		if (sqrSum > 0)
			return sqrt(sqrSum);
		else
			return 0;
	}

	static int elementWisePower(matrix A, double power, matrix *Apower)
	{
		Apower->m = A.m;
		Apower->n = A.n;
		NEW_MAT((*Apower));

		int i;

		for (i = 0; i<A.n*A.m; i++)
			Apower->data[i] = pow(A.data[i], power);

		return EXIT_SUCCESS;
	}

	static int matrixAppendColumns(matrix A, matrix B, matrix *AB)
	{
		if (A.m != B.m)
		{

#ifdef OFFLINE
			printf("Error: append matrixes column wise: the number of row nust be the same\n");
#endif

			return EXIT_FAIL;
		}

		AB->m = A.m;
		AB->n = A.n + B.n;
		NEW_MAT((*AB));

		int i, j;

		for (i = 0; i<A.m; i++)
		{
			for (j = 0; j<A.n + B.n; j++)
			{
				if (j<A.n)
				{
					EL((*AB), i, j) = EL(A, i, j);
				}
				else{
					EL((*AB), i, j) = EL(B, i, j - A.n);
				}
			}
		}

		return EXIT_SUCCESS;
	}

	static int matrixCopy(matrix A, matrix *B)
	{
		B->m = A.m;
		B->n = A.n;
		NEW_MAT((*B));

		int i, j;

		for (i = 0; i<A.m; i++)
		{
			for (j = 0; j<A.n; j++)
			{
				EL((*B), i, j) = EL(A, i, j);
			}
		}

		return EXIT_SUCCESS;
	}

	static int matrixGetColumn(matrix A, int column, matrix *colMat)
	{
		if (column > A.n)
		{
#ifdef OFFLINE
			printf("Error: get a column from a matrix; column to get is larger then the number of columns in the matrix!\n");
#endif

			return EXIT_FAIL;
		}

		colMat->m = A.m;
		colMat->n = 1;
		NEW_MAT((*colMat));

		int i;

		for (i = 0; i<A.m; i++)
		{
			EL((*colMat), i, 0) = EL(A, i, column);
		}

		return EXIT_SUCCESS;
	}

	static int matrixGetRow(matrix A, int row, matrix *rowMat)
	{
		if (row > A.m)
		{
#ifdef OFFLINE
			printf("Error: get a row from a matrix; row index to get is larger then the number of row in the matrix!\n");
#endif
			return EXIT_FAIL;
		}

		rowMat->m = 1;
		rowMat->n = A.n;
		NEW_MAT((*rowMat));

		int i;

		for (i = 0; i<A.n; i++)
		{
			EL((*rowMat), 0, i) = EL(A, row, i);
		}

		return EXIT_SUCCESS;
	}

	static int appendPowerColumnsToMat(matrix A, int maxPower, int beginWithPowerZero, matrix *Apower)
	{

		// Case: beginWithPowerZero = false
		// if A=[x1,x2] and maxPower = 3, this will return [x1,x2,x1^2,x2^2,x1^3,x2^3]
		// Case: beginWithPowerZero = true
		// if A=[x1,x2] and maxPower = 3, this will return [1,1,x1,x2,x1^2,x2^2,x1^3,x2^3]

		if (beginWithPowerZero == false)
		{
			Apower->n = A.n*maxPower;
		}
		else{
			Apower->n = A.n*maxPower + 1;
		}

		Apower->m = A.m;
		NEW_MAT((*Apower));

		int i, j, power;

		for (i = 0; i<A.m; i++)
		{
			for (j = 0; j<A.n; j++) // and NOT (j=0;j<A.n*maxPower;j++)
			{
				if (beginWithPowerZero == false)
				{
					for (power = 0; power<maxPower; power++)
					{
						EL((*Apower), i, j + power*A.n) = pow(EL(A, i, j), power + 1);
					}
				}
				else{
					for (power = 0; power <= maxPower; power++)
					{
						EL((*Apower), i, j + power*A.n) = pow(EL(A, i, j), power);
					}
				}
			}
		}


		return EXIT_SUCCESS;

	}

	static int matrixDifference(matrix A, matrix B, matrix *D)
	{
		if ((A.m != B.m) || A.n != B.n)
		{
#ifdef OFFLINE
			printf("Error: matrix difference: matrix dimensions must be the same");
#endif
			return EXIT_FAIL;
		}

		D->m = A.m;
		D->n = A.n;
		NEW_MAT((*D));

		double v1, v2;
		int i;
		int j;

		for (i = 0; i<A.m; i++)
		{
			for (j = 0; j<A.n; j++)
			{
				v1 = EL(A, i, j);
				v2 = EL(B, i, j);
				EL((*D), i, j) = v1 - v2;
			}
		}

		return EXIT_SUCCESS;
	}

	static int matrixMultiplication(matrix A, matrix B, matrix *AtimesB)
	{
		if (A.n != B.m)
		{
#ifdef OFFLINE
			printf("Error: matrix multiplication: the number of columns of the left matrix must be the same as the number of row in the right matrix\n");
#endif
			return EXIT_FAIL;
		}

		AtimesB->m = A.m;
		AtimesB->n = B.n;
		NEW_MAT((*AtimesB));

		int i, j, k;
		double sum;

		for (i = 0; i<A.m; i++) // a row of A
		{
			for (j = 0; j<B.n; j++) // a column of B
			{
				sum = 0;
				for (k = 0; k<A.n; k++)
				{
					sum += EL(A, i, k) * EL(B, k, j);
				}
				EL((*AtimesB), i, j) = sum;
			}
		}

		return EXIT_SUCCESS;
	}

	static int choleskyDecomposition(matrix X, matrix *L) // A MUST be a square matrix of size shape-by-shape

	{
		if (X.m != X.n)
		{
			return EXIT_FAIL;
		}

		int i, j, k; // loop variables

		L->m = X.m;
		L->n = X.n;
		NEW_MAT((*L));
		double tmp;

		double sum = 0.0;

		for (i = 0; i<X.n; i++)
		{
			for (k = 0; k<i + 1; k++)
				sum += pow(EL((*L), i, k), 2);

			tmp = EL(X, i, i) - sum;

			if (tmp < 0)
			{
				return EXIT_FAIL;
			}

			if (tmp < 0) tmp = 0;
			EL((*L), i, i) = sqrt(tmp);

			sum = 0.0;

			for (j = i + 1; j<X.n; j++)
			{
				for (k = 0; k<i + 1; k++)
					sum += EL((*L), i, k) * EL((*L), j, k);

				EL((*L), j, i) = (EL(X, j, i) - sum) / (EL((*L), i, i) + 1e-30);

				sum = 0.0;
			}
		}

		return EXIT_SUCCESS;
	}

	static int matrixTranspose(matrix A, matrix *AT)
	{
		AT->m = A.n;
		AT->n = A.m;
		NEW_MAT((*AT));

		int i, j;

		for (i = 0; i<A.m; i++)
		{
			for (j = 0; j<A.n; j++)
			{
				EL((*AT), j, i) = EL(A, i, j);
			}
		}

		return EXIT_SUCCESS;
	}

	static int ATransposeTimesB(matrix A, matrix B, matrix *ATB)
	{
		if (A.m != B.m)
		{
			return EXIT_FAIL;
		}

		ATB->m = A.n;
		ATB->n = B.n;
		NEW_MAT((*ATB));

		int i, j, k;
		double sum = 0;

		for (i = 0; i<A.n; i++)
		{
			for (j = 0; j<B.n; j++)
			{
				sum = 0;
				for (k = 0; k<A.m; k++)
					sum += EL(A, k, i) * EL(B, k, j); // we multiply column i with column j

				EL((*ATB), i, j) = sum;

			}
		}

		return EXIT_SUCCESS;
	}

	static int ATimesBTranspose(matrix A, matrix B, matrix *ABT)
	{
		if (A.n != B.n)
		{
			return EXIT_FAIL;
		}

		ABT->m = A.m;
		ABT->n = B.m;
		NEW_MAT((*ABT));

		int i, j, k;
		double sum = 0;

		for (i = 0; i<A.m; i++)
		{
			for (j = 0; j<B.m; j++)
			{
				sum = 0;
				for (k = 0; k<A.n; k++)
					sum += EL(A, i, k) * EL(B, j, k); // we multiply column i with column j

				EL((*ABT), i, j) = sum;

			}
		}

		return EXIT_SUCCESS;
	}

	static int performBackSubstitution(matrix L, matrix u, matrix *w) // given a upper triangle matrix [L^T], solve [L^T][w] = [u]
		// Since our implementation of the cholesky decomposition
		// only calculates the lower triangle matrix [L], we will
		// perform the calculation here by simply reversing the
		// indexes in [L]. THUS HERE [L] IS A LOWER TRIANGLE MATRIX
	{
		if (L.m != u.m)
		{
			return EXIT_FAIL;
		}

		w->m = L.n;
		w->n = u.n;
		NEW_MAT((*w));

		int i, j, k;
		double sum;

		for (k = 0; k<u.n; k++)
		{
			for (i = L.m - 1; i >= 0; i--)
			{
				sum = 0.0;
				for (j = L.n - 1; i<j; j--)
				{
					sum += EL(L, j, i) * EL((*w), j, k);
				}
				EL((*w), i, k) = (EL(u, i, k) - sum) / (EL(L, i, i) + 1e-30);
			}
		}

		return EXIT_SUCCESS;
	}

	static int performForwardSubstitution(matrix L, matrix u, matrix *w) // given a lower triangle matrix [L], solve [L][w] = [u] 
	{
		if (L.m != u.m)
		{
			return EXIT_FAIL;
		}

		w->m = L.n;
		w->n = u.n;
		NEW_MAT((*w));

		int i, j, k;
		double sum;

		for (k = 0; k<u.n; k++)
		{
			for (i = 0; i<L.m; i++)
			{
				sum = 0.0;
				for (j = 0; j<i; j++)
				{
					sum += EL(L, i, j) * EL((*w), j, k);
				}
				EL((*w), i, k) = (EL(u, i, k) - sum) / (EL(L, i, i) + 1e-30);
			}
		}

		return EXIT_SUCCESS;
	}

	static int linearLeastSquares(matrix X, matrix Y, matrix *b)
	{
		if (X.n > X.m)
		{
			return EXIT_FAIL;
		}

		if (X.m != Y.m)
		{
			return EXIT_FAIL;
		}

		// lets first compute [X^T][X] and [X^T][Y]
		matrix XTtimesX, XTtimesY;
		EVAL(ATransposeTimesB(X, X, &XTtimesX));
		EVAL(ATransposeTimesB(X, Y, &XTtimesY));

		// lets then calculate the cholesky decomposition of [X^T][X]
		matrix L;
		EVAL(choleskyDecomposition(XTtimesX, &L));

		matrix w;
		EVAL(performForwardSubstitution(L, XTtimesY, &w));
		EVAL(performBackSubstitution(L, w, b));

		free(XTtimesX.data);
		free(XTtimesY.data);
		free(L.data);
		free(w.data);

		return EXIT_SUCCESS;

	}

	static int polyFit(matrix x, matrix y, int order, int throughZero, matrix *params)
	{
		if (x.n != 1 || y.n != 1)
		{
			return EXIT_FAIL;
		}

		int appendZeroOrder;
		if (throughZero)
			appendZeroOrder = false;
		else
			appendZeroOrder = true;

		matrix xPowers;
		EVAL(appendPowerColumnsToMat(x, order, appendZeroOrder, &xPowers));
		EVAL(linearLeastSquares(xPowers, y, params));

		free(xPowers.data);

		return EXIT_SUCCESS;

	}

	static int computeJacobian(matrix x, matrix params, double(*derivatives[])(matrix, matrix), matrix *jacobian)
	{
		int M = x.m;

		if (jacobian->m != M)
		{
			return EXIT_FAIL;
		}

		if (jacobian->n != params.m)
		{
			return EXIT_FAIL;
		}

		int i, j; // loop veriables
		matrix xi;
		double d;

		for (i = 0; i<M; i++)
		{
			EVAL(matrixGetRow(x, i, &xi));

			for (j = 0; j<params.m; j++)
			{
				d = (derivatives[j])(xi, params);
				EL((*jacobian), i, j) = d;
			}
			free(xi.data);
		}

		return EXIT_SUCCESS;
	}

	static int iteration(matrix x,
		matrix y,
		matrix *params,
		int(*fitFunction)(matrix, matrix, matrix*),
		double(*derivatives[])(matrix, matrix),
		double *precesion)
	{

		int N = x.m;
		int nParams = params->m;
		matrix fit;
		fit.m = N;
		fit.n = 1;
		NEW_MAT(fit);

		fitFunction(x, *params, &fit);
		matrix r;
		EVAL(matrixDifference(y, fit, &r));


		matrix jacobian;
		jacobian.m = N;
		jacobian.n = nParams;
		NEW_MAT(jacobian);
		EVAL(computeJacobian(x, *params, derivatives, &jacobian))

			// now solve jacobian * u = r
			matrix u;
		EVAL(linearLeastSquares(jacobian, r, &u));

		*precesion = vectorMagnitude(u);

		int i;
		for (i = 0; i<params->m; i++){ params->data[i] -= u.data[i]; }

		// free all the meory used
		free(fit.data);
		free(r.data);
		free(jacobian.data);
		free(u.data);

		return EXIT_SUCCESS;
	}

	static int gaussNewton(matrix x,
		matrix y,
		matrix paramsInit,
		int(*fitFunction)(matrix, matrix, matrix*),
		double(*derivatives[])(matrix, matrix),
		matrix *paramsFinal)
	{
		double precision = 0.001;
		int countIter = 0, maxIter = 100;
		double d = 1.0;
		EVAL(matrixCopy(paramsInit, paramsFinal));

		int exitCode = EXIT_SUCCESS;

		while (d > precision && exitCode == EXIT_SUCCESS)
		{
			countIter += 1;
			if (iteration(x, y, paramsFinal, fitFunction, derivatives, &d) == EXIT_FAIL)
			{
				exitCode = EXIT_FAIL;
			}

			if (countIter == maxIter)
			{
				exitCode = EXIT_FAIL;
			}
		}

		return exitCode;
	}


	static int fitFunction(matrix xy, matrix params, matrix *modelValues)
	{
		int M = xy.m;
		int i;

		double mux = EL(params, 0, 0);
		double muy = EL(params, 1, 0);

		double sigx = EL(params, 2, 0);
		double sigy = EL(params, 3, 0);

		double A = EL(params, 4, 0);
		double z0 = EL(params, 5, 0);
		double xi, yi, zi;

		double wx, wy;

		for (i = 0; i<M; i++)
		{
			xi = EL(xy, i, 0);
			yi = EL(xy, i, 1);

			wx = (xi - mux) / ((sqrt(2.0f)*sigx) + 1e-30);
			wy = (yi - muy) / ((sqrt(2.0f)*sigy) + 1e-30);

			zi = A*exp(-pow(wx, 2) - pow(wy, 2)) + z0;

			EL((*modelValues), i, 0) = zi;
		}

		return EXIT_SUCCESS;
	}

	static double derivative_ROverMux(matrix xy, matrix params)
	{
		double xi = EL(xy, 0, 0), yi = EL(xy, 0, 1);

		double mux = EL(params, 0, 0);
		double muy = EL(params, 1, 0);

		double sigx = EL(params, 2, 0);
		double sigy = EL(params, 3, 0);

		double A = EL(params, 4, 0);

		double wx = (xi - mux) / ((sqrt(2.0)*sigx) + 1e-30);
		double wy = (yi - muy) / ((sqrt(2.0)*sigy) + 1e-30);

		return -2 * wx*A*exp(-pow(wx, 2) - pow(wy, 2)) / ((sqrt(2.0)*sigx) + 1e-30);
	}

	static double derivative_ROverMuy(matrix xy, matrix params)
	{
		double xi = EL(xy, 0, 0), yi = EL(xy, 0, 1);

		double mux = EL(params, 0, 0);
		double muy = EL(params, 1, 0);

		double sigx = EL(params, 2, 0);
		double sigy = EL(params, 3, 0);

		double A = EL(params, 4, 0);

		double wx = (xi - mux) / ((sqrt(2.0)*sigx) + 1e-30);
		double wy = (yi - muy) / ((sqrt(2.0)*sigy) + 1e-30);


		return -2 * wy*A*exp(-pow(wx, 2) - pow(wy, 2)) / ((sqrt(2.0)*sigy) + 1e-30);
	}

	static double derivative_ROverSigx(matrix xy, matrix params)
	{
		double xi = EL(xy, 0, 0), yi = EL(xy, 0, 1);
		double mux = EL(params, 0, 0);
		double muy = EL(params, 1, 0);
		double sigx = EL(params, 2, 0);
		double sigy = EL(params, 3, 0);

		double A = EL(params, 4, 0);

		double wx = (xi - mux) / ((sqrt(2.0)*sigx) + 1e-30);
		double wy = (yi - muy) / ((sqrt(2.0)*sigy) + 1e-30);

		double tmp = (sqrt(2.0)*pow(sigx, 2))*exp(-pow(wx, 2) - pow(wy, 2));
		if (tmp == 0) return 0;

		return -2 * A*wx*(xi - mux) / tmp;
	}

	static double derivative_ROverSigy(matrix xy, matrix params)
	{
		double xi = EL(xy, 0, 0), yi = EL(xy, 0, 1);

		double mux = EL(params, 0, 0);
		double muy = EL(params, 1, 0);

		double sigx = EL(params, 2, 0);
		double sigy = EL(params, 3, 0);

		double A = EL(params, 4, 0);

		double wx = (xi - mux) / ((sqrt(2.0)*sigx) + 1e-30);
		double wy = (yi - muy) / ((sqrt(2.0)*sigy) + 1e-30);

		double tmp = (sqrt(2.0)*pow(sigy, 2))*exp(-pow(wx, 2) - pow(wy, 2));
		if (tmp == 0) return 0;

		return -2 * A*wy*(yi - muy) / tmp;
	}

	static double derivative_ROverA(matrix xy, matrix params)
	{
		// differentiate 
		// z = A*exp(-wx^2 - wy^2) + z0, with wx = (x-mux)/(sqrt(2)*sigx), wy = (y-muy)/(sqrt(2)*sigy)
		// with respect to A (we return the negative derivative -dz/dA !)

		double xi = EL(xy, 0, 0), yi = EL(xy, 0, 1);

		double mux = EL(params, 0, 0);
		double muy = EL(params, 1, 0);

		double sigx = EL(params, 2, 0);
		double sigy = EL(params, 3, 0);

		double wx = (xi - mux) / (sqrt(2.0)*sigx + 1e-30);
		double wy = (yi - muy) / (sqrt(2.0)*sigy + 1e-30);

		return -exp(-pow(wx, 2) - pow(wy, 2));
	}

	static double derivative_ROverZnull(matrix x, matrix params)
	{
		return -1;
	}

	static int estimateParams(matrix xy, matrix z, matrix *estimates)
	{
		/*
		An improvement of this function would be to sort the z-values and take the mean of the bottom x%
		as minZ
		*/

		int i, M = xy.m;
		double mux = 0, muy = 0;
		double sigx = 0, sigy = 0;
		double xi, yi, zi, totalZ = 0;
		double maxY = -1E6, minY = 1E6;
		double maxX = -1E6, minX = 1E6;
		double maxZ = -1E6, minZ = 1E6;

		for (i = 0; i<M; i++)
		{
			xi = EL(xy, i, 0);
			yi = EL(xy, i, 1);
			zi = EL(z, i, 0);

			if (xi > maxX) { maxX = xi; }
			if (xi < minX) { minX = xi; }
			if (yi > maxY) { maxY = yi; }
			if (yi < minY) { minY = yi; }
			if (zi > maxZ) { maxZ = zi; }
			if (zi < minZ) { minZ = zi; }
		}

		double A = maxZ - minZ;

		for (i = 0; i<M; i++)
		{
			xi = EL(xy, i, 0);
			yi = EL(xy, i, 1);
			zi = EL(z, i, 0) - minZ;

			if (zi < 0.2*A)
				continue;

			totalZ += zi;
			mux += xi*zi;
			muy += yi*zi;
		}

		mux /= (totalZ + 1e-30);
		muy /= (totalZ + 1e-30);

		if (M < 0) M = 0;
		double dx = (maxX - minX) / (sqrt((float)M) + 1e-30); // we assume that the sampling distance in x and y is approx.
		double dy = (maxY - minY) / (sqrt((float)M) + 1e-30); // the same, so the number of samples along x and y is ~sqrt(M)

		double volumeUnderGauss = totalZ*dx*dy;

		if (volumeUnderGauss < 0) volumeUnderGauss = 0;
		double tmp = 0;
		if (A != 0) tmp = volumeUnderGauss / (A * 1.0) / (2 * 3.14);
		if (tmp < 0) tmp = 0;
		sigx = sqrt(tmp);
		sigy = sigx;

		estimates->m = 6;
		estimates->n = 1;
		NEW_MAT((*estimates));
		estimates->data[0] = mux;
		estimates->data[1] = muy;
		estimates->data[2] = sigx;
		estimates->data[3] = sigy;
		estimates->data[4] = A;
		estimates->data[5] = minZ;

		return EXIT_SUCCESS;
	}

	static int gaussian2DFitn(double x[],
		double y[],
		double intensities[],
		int dataSize,
		double* paramsArray)
	{
		int result = EXIT_SUCCESS;

		matrix xy, z;
		xy.m = dataSize;
		xy.n = 2;
		NEW_MAT(xy);

		z.m = dataSize;
		z.n = 1;
		NEW_MAT(z);

		int i;
		for (i = 0; i<dataSize; i++)
		{
			EL(xy, i, 0) = x[i];
			EL(xy, i, 1) = y[i];
			EL(z, i, 0) = intensities[i];
		}

		matrix paramsInit, params;
		EVAL(estimateParams(xy, z, &paramsInit));

		double(*derivatives[])(matrix, matrix) = { derivative_ROverMux,
			derivative_ROverMuy,
			derivative_ROverSigx,
			derivative_ROverSigy,
			derivative_ROverA,
			derivative_ROverZnull };

		if (gaussNewton(xy, z, paramsInit, fitFunction, derivatives, &params) == EXIT_FAIL)
		{

			result = EXIT_FAIL;
		}
		else{

			for (i = 0; i < 6; i++){ if (params.data[i] == params.data[i]) paramsArray[i] = params.data[i]; else  paramsArray[i] = 0; }

		}

		free(params.data);
		free(paramsInit.data);
		free(xy.data);
		free(z.data);

		return result;

	}
	}

#endif 